---
title: "GFP32_Slingshot for Oligodendroglia"
author: "Nina-Lydia Kazakou"
date: "19/03/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Set-up

### Load libraries
```{r message=FALSE, warning=TRUE}
library(slingshot)
library(scales)
library(Seurat)
library(ggplot2)
library(RColorBrewer)
library(ggsci)
library(here)
library(tidymodels)
```

### Colour Palette

```{r load_palette}
mypal <- pal_npg("nrc", alpha = 0.7)(10)
mypal2 <-pal_tron("legacy", alpha = 0.7)(7)
mypal3<-pal_lancet("lanonc", alpha = 0.7)(9)
mypal4<-pal_simpsons(palette = c("springfield"), alpha = 0.7)(16)
mypal5 <- pal_rickandmorty(palette = c("schwifty"), alpha = 0.7)(6)
mypal6 <- pal_futurama(palette = c("planetexpress"), alpha = 0.7)(5)
mypal7 <- pal_startrek(palette = c("uniform"), alpha = 0.7)(5)
mycoloursP<- c(mypal, mypal2, mypal3, mypal4)
```


### Load Object 
```{r}
GFP32 <- readRDS(here("data", "Processed", "Original_Objects", "GFP32_annotated_srt.rds"))
Idents(GFP32) <- "Cluster_id"
```

```{r}
DimPlot(GFP32, group.by = "Cluster_id", pt.size = 0.9, label = TRUE, cols = mycoloursP) & NoLegend()
```

```{r Oligodendroglia_Subset}
oligos <- subset(GFP32, idents = c("OPC 1", "OPC 2", "OLIGO 1", "OLIGO 2", "OLIGO 3"))
```

```{r fig.height=5, fig.width=8}
DimPlot(oligos, group.by = "Cluster_id", split.by = "orig.ident", pt.size = 0.9, label = TRUE, cols = mycoloursP) & NoLegend()
```

```{r}
sce <- as.SingleCellExperiment(oligos)
```

```{r eval=FALSE}
dir.create(here("outs", "Processed", "GFP32", "Pseudo-time_Analysis"))
dir.create(here("outs", "Processed", "GFP32", "Pseudo-time_Analysis", "Slingshot"))
```

# Slingshot
Running slingshot is a two-step process composed of identifying the global lineage structure with a cluster-based minimum spanning tree (MST) and fitting simultaneous principal curves to describe each lineage.

```{r}
# Use wrapper function 
sce <- slingshot(sce, clusterLabels = 'Cluster_id', reducedDim = 'UMAP')
```

The output is an sce with slingshot results in the colData. All of the results are stored in a PseudotimeOrdering object, which is added to the colData of the original object and be accessed with colData(sce)$slingshot. 
All the inferred pseudotime variables (one per lineage) are added to the colData, individually. 

To extract all slingshot results in a single object, we will use SlingshotDataSet, objects that are primarily used for visualization, as several plotting methods are included with the package. 
```{r}
summary(sce$slingPseudotime_1)
```

```{r fig.height=6, fig.width=8}
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)

plotcol <- colors[cut(sce$slingPseudotime_1, breaks=100)]

plot(reducedDims(sce)$UMAP, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, col='black')
```
We can also see how the lineage structure was initially estimated by the cluster-based minimum spanning tree:
```{r fig.height=6, fig.width=8}
plot(reducedDims(sce)$UMAP, col = mycoloursP[sce$Cluster_id], pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, type = 'lineages', col = 'black')
```

Slingshot also allows for the specification of known endpoints.
Clusters which are specified as terminal cell states will be constrained to have only one connection when the MST is constructed (i.e., they must be leaf nodes). This constraint could potentially impact how other parts of the tree are drawn.

```{r}
sds <- slingshot(Embeddings(oligos, "umap"), clusterLabels = oligos@meta.data$Cluster_id, start.clus = "OPC 1") # set pri-OPC as the starting cluster for the pseudotime

cl <- oligos@meta.data$Cluster_id
cols <- mycoloursP

plot(slingReducedDim(sds), col = cols[cl], pch = 16, cex = 0.5) 
lines(SlingshotDataSet(sds), lwd = 2, type = 'lineages', col = 'black') 
```

I will also try Slingshot on different time-points, since I've noticed in the past that the trajectory is altered between oligodendroglia that has been cultured for one week and oligodendroglia that has been cultured for longer. 

```{r}
Idents(oligos) <- "orig.ident"
wk_1 <- subset(oligos, ident = "GFP32_wk1")
wk_2 <- subset(oligos, ident = "GFP32_wk2")
wk_3 <- subset(oligos, ident = "GFP32_wk3")
```

```{r}
wk_1_sce <- as.SingleCellExperiment(wk_1)
wk_2_sce <- as.SingleCellExperiment(wk_2)
wk_3_sce <- as.SingleCellExperiment(wk_3)
```

```{r}
wk_1_sce <- slingshot(wk_1_sce, clusterLabels = 'Cluster_id', reducedDim = 'UMAP')

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)

plotcol <- colors[cut(wk_1_sce$slingPseudotime_1, breaks=100)]

plot(reducedDims(wk_1_sce)$UMAP, col = mycoloursP[wk_1_sce$Cluster_id], pch=16, asp = 1)
lines(SlingshotDataSet(wk_1_sce), lwd=2, type = 'lineages', col = 'black')
```


```{r}
wk_2_sce <- slingshot(wk_2_sce, clusterLabels = 'Cluster_id', reducedDim = 'UMAP')

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)

plotcol <- colors[cut(wk_2_sce$slingPseudotime_1, breaks=100)]

plot(reducedDims(wk_2_sce)$UMAP, col = mycoloursP[wk_2_sce$Cluster_id], pch=16, asp = 1)
lines(SlingshotDataSet(wk_2_sce), lwd=2, type = 'lineages', col = 'black')
```

```{r}
wk_3_sce <- slingshot(wk_3_sce, clusterLabels = 'Cluster_id', reducedDim = 'UMAP')

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)

plotcol <- colors[cut(wk_3_sce$slingPseudotime_1, breaks=100)]

plot(reducedDims(wk_3_sce)$UMAP, col = mycoloursP[wk_3_sce$Cluster_id], pch=16, asp = 1)
lines(SlingshotDataSet(wk_3_sce), lwd=2, type = 'lineages', col = 'black')
```


```{r}
sessionInfo()
```